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We  report  on  the  implementation  of  Burgers  equation  as  a  type-II  quantum  computation  on 
an  NMR  quantum  information  processor.  Since  the  flow  field  evolving  under  the  Burgers  equation 
develops  sharp  features  over  time,  this  is  a  better  test  of  liquid  state  NMR  implementations  of  type-II 
quantum  computers  than  the  previous  examples  using  the  diffusion  equation.  In  particular,  we  show 
that  Fourier  approximations  used  in  the  encoding  step  are  not  the  dominant  error.  Small  systematic 
errors  in  the  collision  operator  accumulate  and  swamp  all  other  errors.  We  propose,  and  demonstrate, 
that  the  accumulation  of  this  error  can  be  avoided  to  a  large  extent  by  replacing  the  single  collision 
operator  with  a  set  of  operators  with  random  errors  and  similar  fidelities.  Experiments  have  been 
implemented  on  16  two-qubit  sites  for  eight  successive  time  steps  for  the  Burgers  equation. 

PACS  numbers: 


It  has  been  suggested  that  some  classical  computa¬ 
tional  problems  can  be  solved  by  using  a  hybrid  clas¬ 
sical  quantum  device,  a  type  II  quantum  computer  [1,  2]. 
Such  a  device  is  essentially  an  array  of  small  quantum  im- 
formation  processors  (QIP)  sharing  information  through 
classical  channels.  NMR  has  proven  to  be  a  useful  testbed 
for  QIP,  and  in  particular  we  have  shown  that  a  lattice 
of  parallel  QIPs  can  be  mapped  onto  a  spin  system  by 
creating  a  correspondence  between  the  lattice  sites  and 
spatially  distinct  spin  ensembles.  A  first  proof-of-concept 
for  numerically  predicting  the  time-dependent  solution 
of  classical  partial  differential  equation  with  dissipative 
terms  using  our  NMR  technique  was  demonstrated  for 
the  diffusion  equation  [3,  4]. 

One  of  the  most  important  challenges  to  implement¬ 
ing  a  useful  type-II  quantum  architecture  is  to  avoid  the 
accumulation  of  systematic  errors.  In  the  NMR  imple¬ 
mentations  to  date  there  are  two  important  sources  of 
systematic  errors:  (1)  a  linear  approximation  relating  the 
excited  magnetization  to  the  Fourier  components  of  the 
shaped  RF  pulse;  and  (2)  errors  from  the  repeated  col¬ 
lision  operators.  Here  we  explore  the  impact  of  these 
errors  on  a  simple  computation  and  illustrate  a  simple 
means  of  reducing  the  accumulated  error. 

The  ensemble  nature  of  the  spin  system  allows  us 
to  split  the  sample  into  a  spatial  array  of  lattice  sites. 
Well  developed,  methods  from  magnetic  resonance  imag¬ 
ing  (MRI)  allow  us  to  selectively  address  the  spins  in 
each  of  these  sites.  Typically  the  addressing  is  carried 
out  in  a  space  reciprocal  to  the  spatial  mapping,  called 
fc-space,  where  fc  is  the  wave-number  of  the  correspond¬ 
ing  Fourier  components.  The  fc-space  formalism  [5]  pro¬ 
vides  a  recipe  for  writing  a  spatially  varying  spin  rotation 
across  an  ensemble  of  spins  that  have  been  distinguished 
from  each  other  by  a  magnetic  field  gradient.  The  fc- 
space  formalism  is  essentially  the  application  of  shaped 
radio  frequency  (RF)  pulses  in  the  presence  of  a  linear 
magnetic  gradient  field  as  a  means  of  exciting  selective 
frequencies.  For  most  studies  the  full  k-space  formalism 
is  not  employed  and  a  linear  approximation  is  invoked.  If 


the  rotation  angle  of  the  shaped  pulse  is  small,  then  the 
excited  magnetization  may  be  accurately  calculated  only 
to  first  order  in  that  angle,  and  the  excited  magnetiza¬ 
tion  is  related  to  the  RF  waveform  simply  by  a  Fourier 
transform.  As  a  result,  the  required  RF  waveform  can 
also  be  determined  by  taking  the  inverse  Fourier  trans¬ 
form  of  the  desired  initial  magnetization.  This  technique 
allows  us  to  encode  arbitrary  magnetization  profiles  span¬ 
ning  the  various  spatial  locations  in  our  experiment  and 
thereby  approximating  any  desired  initial  conditions.  In 
the  previously  implemented  diffusion  equation,  higher  or¬ 
der  Fourier  components  of  the  number  density  are  atten¬ 
uated  by  the  dynamics  and  the  solution  is  stable  even  in 
the  presence  of  substantial  accumulated  errors. 

To  push  the  development  of  type-II  implementations 
we  have  chosen  to  explore  the  nonlinear  Burgers  equa¬ 
tion  to  test  the  breakdown  for  the  Fourier  approximation. 
Over  time,  a  shock  front  forms  and  high  spatial  frequen¬ 
cies  in  the  magnetization  profile  become  important  and 
it  is  these  high  spatial  frequencies  that  we  expect  to  be 
most  sensitive  to  errors.  The  numerical  treatment  of  the 
QLG  algorithm  for  the  Burgers  equation  [6]  therefore  of¬ 
fers  a  stronger  proof  of  our  NMR  quantum  computing  ap¬ 
proach  since  the  effect  of  the  nonlinear  convective  term  in 
the  equation  generates  a  sharp  edge  as  a  shock  develops 
in  time  that  is  not  mimicked  by  spin  relaxation,  random 
self-diffusion,  nor  RF  inhomogeneities.  In  addition,  we 
demonstrate  shock-formation  driven  by  a  tunable  viscos¬ 
ity  parameter  to  show  that  the  width  of  the  shock  front 
is  not  determined  by  implementation  imperfections. 

The  first-order  accurate  Fourier  approximation  was  ex¬ 
pected  to  be  the  dominant  error  source  in  the  NMR 
implementation.  However,  NMR  simulations  with  con¬ 
trolled  errors  shows  that  the  systematic  error  induced  by 
the  experimental  implementation  of  the  unitary  collision 
operator  associated  with  the  quantum  lattice  gas  (QLG) 
algorithm  is  the  major  challenge.  Replacing  the  single 
collision  operator  with  a  set  of  operators  to  randomize 
errors  allows  us  to  improve  the  robustness  of  the  imple¬ 
mentation. 


20050829  025 


2 


.i. 


A  quantum  lattice  gas  is  a  system  of  quantum  particles 
moving  and  colliding  on  a  discrete  spacetime  lattice.  This 
quantum  particle  system  is  isomorphic  to  a  lattice-based 
qubit  system.  The  mapping  is  as  follows:  the  probabil¬ 
ity  of  a  particle  residing  at  a  particular  lattice  node  is 
equated  to  the  moduli  squared  of  the  probability  ampli¬ 
tude  of  a  qubit  at  a  unique  location  being  in  its  excited 
quantum  state.  That  is,  each  spatial  location  that  a  par¬ 
ticle  may  occupy  is  mapped  onto  a  qubit  associated  with 
a  unique  location. 

The  dynamics  of  evolution  in  the  QLG  algorithm  can 
be  described  in  three  scales,  the  microscopic,  mesoscopic, 
and  macroscopic  scales.  At  the  microscopic  scale,  each 
particle  has  some  probability  of  moving  along  the  lat¬ 
tice.  For  example,  the  particle  can  move  to  the  right 
or  left  lattice  site  in  a  one-dimentional  construction.  A 
simplified  dynamics  allows  a  particle  to  change  its  direc¬ 
tion  of  motion  (via  a  collision  with  another  particle)  or 
keep  moving  at  constant  speed  in  its  original  direction 
of  motion  ( streaming ).  One  can  describe  the  dynamical 
behavior  of  particles  at  the  mesosopic  scale  by  determin¬ 
ing  their  occupation  probablities  on  the  lattice  points;  in 
one-dimention,  only  left-moving  and  right-moving  prob¬ 
abilities  are  needed.  The  mesoscopic  dynamical  behav¬ 
ior  of  the  system  is  modeled  by  a  finite-difference  form 
of  a  quantum  Boltzmann  equation.  Finally,  to  bridge 
to  the  macroscopic  scale,  the  occupation  probabilities  of 
the  particles  residing  at  each  lattice  site  are  summed  to¬ 
gether  to  determine  the  number  density.  This  number 
density  quantity  defined  at  each  lattice  node  becomes  a 
continuous  field  at  the  lattice  resolution  approaches  in¬ 
finity,  which  is  called  the  continuum  limit.  Through  a 
Chapman-Enskog  purturbation  procedure  applied  in  the 
continuum  limit,  from  the  quantum  Boltzmann  equation 
emerges  an  effective  field  theory  that  is  parabolic  in  time 
and  space  and  nonlinear  in  the  number  density  [7]. 

The  QLG  algorithm  is  initialized,  in  the  NMR  case,  by 
encoding  the  particles’  occupation  probabilities  as  a  spin- 
magnetization  profile.  To  handle  the  one-dimensional 
Burgers  equation,  it  is  sufficient  to  use  two  qubits  (two 
spin- 1  nuclei)  per  lattice  site,  where  each  stores  a  single 
real  valued  occupation  probability.  A  room-temperature 
solution  of  isotopically-labeled  chloroform  ( 13CHClz ) 
was  chosen  for  implementing  the  experiments,  where  the 
hydrogen  and  the  labeled  carbon  nucleus  served  as  qubits 
1  and  2,  respectively.  The  difference  of  the  gyro-magnetic 
ratio  of  two  spins  generates  widely  spaced  resonant  fre¬ 
quencies  that  allows  us  to  address  each  spin  indepen¬ 
dently.  A  lattice  of  QIPs  are  related  to  the  ensemble 
sample  by  creating  a  correspondence  between  lattice  sites 
and  spatially  dependent  positions  in  the  sample.  A  lin¬ 
ear  magnetic  field  gradient  is  used  to  generate  distinct 
spatially-dependent  resonant  frequencies  that  we  can  dis¬ 
tinguish  and  modulate  by  a  shaped  RF  pulse.  In  this  way, 
the  magnetic  field  gradient  allows  the  entire  spin  ensem¬ 
ble  to  be  sliced  into  a  lattice  of  smaller,  and  individually 
addressable,  sub-ensembles. 

The  lattice  initialization  starts  by  transforming  ther- 
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FIG.  1:  QLG  algorithm  implemented  in  four  steps.  Three  hori¬ 
zontal  lines  represent  proton  spin,  carbon  spin  and  fried  gradients. 
Both  starting  magnetizations  are  encoded  in  proton  channel  first 
due  to  the  high  signal  to  noise  ratio  while  decoupled  in  carbon  chan¬ 
nel  to  prevent  interfering  of  scalar  coupling.  The  collision  operator 
is  applied  after  the  initialization.  Measurements  are  also  taken  in 
two  steps  in  the  proton  channl  followed  by  data  processing  in  a 
personal  computer. 


mal  equilibrium  states  into  pseudo-pure  states  [8].  The 
equilibrium  state  is  highly  mixed  and  the  two  nuclear 
spins  have  unequal  magnetizations.  Thus,  equalization 
of  the  magnetizations  is  required  prior  to  creating  the 
pseudo-pure  state.  The  dynamical  evolution  is  caused 
by  a  collision  operator  (a  quantum  operation),  and  mea¬ 
surement  and  streaming  (classical  operations)  according 
to  the  QLG  algorithmic  paradigm.  The  four  main  sec¬ 
tions  of  the  NMR  implementation  of  QLG  algorithm  are 
graphically  depicted  in  Figure  1. 

First,  each  occupation  probablility  is  mapped  onto  a 
lattice  site  as  the  expectation  value  of  a  number  opera¬ 
tor  at  a  spacetime  site  nAx  at  time  mAt.  As  a  result, 
the  initial  state  of  the  ath  qubit  is  y/  fa(nAx,mAt)\l)  + 
y/i  —  fa(nAx,mAt)\0).  The  combined  the  wave  func¬ 
tion  for  a  lattice  site  is  a  tensor  product  over  the  qubits, 

\ip(nAx,mAt)}  =  y/Kh\U)  +  i//i(  1  -  /2)|10>  (1) 

+  + y/{i  -  h)(  i  -  /2)ioo> 


In  the  basis  of  a  two-qubit  system,  the  number  opera¬ 
tors  for  the  occupancy  of  qubits  are  defined  in  terms  of 


the  singleton  qubit  number  operation  h  =  ^  ^  ^  J  as 

follows:  hi  =  1  ®  n  and  h2  —  n®  1.  Therefore,  the 
occupation  probability  is  represented  as  follows: 


fa {nAx, mAt)  =  (ip(nAx,mAt)\ha\ip(nAx,mAt)).  (2) 

The  macroscopic  scale  dynamical  quantity  of  the  quan¬ 
tum  lattice  gas  is  the  number  density,  p,  defined  as  the 
sum  of  the  occupancy  probablity.  The  equilibrium  occu¬ 
pation  probabilities  that  we  use  are 


where  ea  is  ±1  for  different  qubits. 

The  initial  magnetization  is  specified  by  using  a  RF 
pulse  shaped  by  the  Fourier  transform  of  the  desired 
magnetization  (tranform  of  the  initial  number  density 
profile).  While  applying  the  shaped  pulse,  a  carbon  de¬ 
coupling  sequence  is  performed  to  prevent  the  scalar  cou¬ 
pling  from  interfering  with  the  low  power  shaped  pulses. 
In  addition,  the  J  pulse,  which  rotates  the  information 
from  the  x-axis  to  the  z-axis,  is  applied  separately  just 
after  each  initialization.  This  is  done  to  keep  the  valu¬ 
able  information  along  the  longitudinal  direction  where 
it  will  not  be  affected  by  the  gradient  and  chemical  shift. 
The  encoding  of  initial  states  on  both  spins  is  accom¬ 
plished  in  two  steps:  The  initial  carbon  magnetication  is 
recorded  on  the  protons  before  being  transferred  to  the 
carbons  and  followed  by  the  initialization  of  proton  mag¬ 
netication.  Furthermore,  a  short  pulse  sequence,  called 
the  clean  sequence,  is  executed  after  the  first  swap  gate 
to  erase  the  phase  distortion  that  may  be  caused  by  the 
decoupling  sequence. 

Second,  the  evolution  of  fa  is  governed  by  the  com¬ 
bined  action  of  the  collision  operator,  measurement  and 
streaming.  The  collision  operator  is  applied  to  all  the 
lattice  sites  independently,  resulting  in  |i/)  ( nAx ))  = 
C\ij)(nAx)),  for  all  n.  The  choice  of  the  particular  com¬ 
ponents  of  the  unitary  collision  operator  determines  the 
form  of  the  macroscopic  effective  field  theory  (a  parabolic 
partial  differential  equation)  and  the  value  of  its  trans¬ 
port  coefficients  (coefficients  of  the  dissipative  terms).  A 
general  representation  of  the  collision  operator  for  the 
Burgers  equation  is  a  block  diagonal  matrix.  This  single 
quantum  operator  is  chosen  to  be 

C  =  exp  ~  tf°x)\ .  (4) 


which  has  the  following  matrix  representation: 

/I  0  0  0\ 

A  j  .0  0.8  0.6  0  I 

0  -0.6  0.8  0  I  - 

\0  0  01/ 


(5) 


The  unitary  operator  C  can  be  decomposed  of  a  se¬ 
quence  of  RF  pulses  and  scalar  coupling.  The  product 
operators  in  the  exponent  commute  with  each  other,  re¬ 
sulting  in  C  =  exp  HjJajO'fof]  «P  Hzimrf'7*}- 
Both  terms  can  be  expanded  as  natural  scalar  Hamilto¬ 
nian  couplings  sandwiched  with  the  appropriate  single 
rotations,  resulting  in 


C  =  e~iHfrv  +(Ty)e-i*kaz  e-*4is2  (6) 


The  exponential  terms  of  single  spin  rotations  are  imple¬ 
mented  by  7t/2  and  7r/4  pulses.  The  exponents  of  terms 
with  crfcrf  represent  the  natural  internal  Hamiltonian 
evolutions  with  time  period  1/2J,  where  J  is  214j Hz. 
Here,  the  evolution  of  the  internal  Hamiltonian  is  ignored 


while  the  RF  pulse  is  applied.  This  approximation  leads 
to  a  systematic  error  that  will  accumulate  during  the 
course  of  the  computation.  In  general,  these  errors  are 
easy  to  avoid,  but  since  the  purpose  of  the  investigation 
was  to  explore  the  sensitivity  to  accumulated  errors  we 
did  not  correct  it.  The  collision  operator  follows  the  en¬ 
coding  (Step  2),  and  it  is  implemented  without  magnetic 
field  gradients  to  ensure  that  all  of  the  sites  in  the  sample 
undergo  the  same  transformation. 

Third,  we  measure  the  occupation  probabilities.  This 
process  erases  all  the  superpositions  and  quantum  entan¬ 
glement  that  was  created  by  the  unitary  collision  opera¬ 
tor  in  the  second  step. 

The  occupation  numbers  of  each  spin  are  ob¬ 
tained  following  the  collision  step  by  measuring  the  z- 
magnetization  according  to  the  following  equation 

fa(n,m)  =  i  [1  +  (ip(n,  m)|cr“|^>(n,  to))].  (7) 

Since  only  ax  and  ay  are  observable  in  our  NMR  spec¬ 
trometer,  a  7r/2  pulse  has  been  used  to  bring  the  z- 
magnetization  into  the  transverse  plane.  The  measure¬ 
ments  are  done  in  two  seperate  experiments,  where  a 
SWAP  gate  is  applied  to  bring  the  magnetization  from 
carbon  channel  to  the  proton  channel.  This  SWAP  op¬ 
eration  is  done  because  the  higher  signal-to-noise  ratio 
in  the  proton  channel  allows  us  to  improve  the  accuracy 
of  our  implementation.  During  the  “readout”  process 
(Step  3),  a  week  magnetic  field  gradient  is  applied  to  dis¬ 
tinguish  different  sites.  The  observed  proton  signals  are 
digitized  and  Fourier  transformed,  allowing  us  to  record 
the  spatially-dependent  spin  magnetization  profile. 

Fourth,  arid  last  step  of  the  QLG  algorithm,  we  shift 
the  fa  obtained  in  the  previous  step  to  its  nearest  neigh¬ 
bor  using  a  short  Matlab  program.  This  step  requires 
only  classical  communication  between  neighboring  sites. 
The  time  is  incremented  after  this  step.  Then,  we  loop 
back  to  step  1  and  update  the  field  of  occupation  proba¬ 
bilities  over  the  lattice  sites.  In  this  way,  we  can  continue 
to  iterate  forward  in  time  and  make  a  time-history  record 
of  the  occupation  probabilities,  which  in  turn  gives  us  the 
temporal  evolution  of  the  number  density  field.  In  the 
implementation  of  the  Burger  equation,  we  observed  de¬ 
viations  between  the  numerically  predicted  data  points 
and  analytically  predicted  solutions.  These  errors  can 
be  attributed  to  imperfections  in  the  NMR  implementa¬ 
tion.  The  major  error  sources  in  the  NMR  implemen¬ 
tation  are  known,  so  to  explore  the  source  and  relative 
strength  of  these  errors,  we  have  simulated  the  NMR  ex¬ 
periments.  The  major  error  source  in  this  implemen¬ 
tation  is  the  collision  operator,  and  it  is  introduced  by 
ignoring  the  scalar  coupling  between  proton  arid  carbon 
during  the  RF  pulses.  When  applying  an  RF  pulse  on 
the  proton  qubit,  the  Hamiltonian  in  the  rotating  form  is 
H  =  2-nJu^a^  +  ,  where  B\  is  the  strength  of 

the  RF  pulse.  With  the  presence  of  the  scalar  coupling, 
a  small  portion  of  the  proton  magnetization  has  been 
transfered  to  the  carbon  qubit.  Therefore,  the  applied 
propagator  can  be  recast  as  U  UdeBiredUetror. 
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FIG.  2:  The  growth  of  the  systematic  errors  due  to  the  collision 
operator  in  two  NMR  implementation.  The  single  collision  opera¬ 
tor  data  (dots)  is  fit  (solid  line)  with  a  line  of  slope  1,  which  shows 
linear  growth  of  the  error.  The  collision  operator  data  with  mod¬ 
ulated  phases  (pluses)  is  the  fit  with  a  line  of  slope  3/4  (dashed 
line).  The  buildup  of  the  systematic  errors  has  been  slowed  down 
by  proposed  method.  However,  the  systematic  errors  have  not  been 
totally  converted  into  random  errors. 


The  error  in  the  collision  operator  is  a  systematic  er¬ 
ror  that  builds  up  throughout  the  successive  time  steps. 
Although  this  is  not  the  dominant  error  at  the  begin¬ 
ning  of  the  implementation,  it  eventually  dominates  the 
first-order  error  due  to  the  Fourier  approximation  and 
becomes  the  dominant  issue  after  just  several  time  step 
interations.  Notice  that  while  the  reduction  of  the  initial 
magnetization  from  the  Fourier  transform  is  systematic, 
since  the  magnetization  profile  is  changing  the  errors  are 
not  precisely  repeated.  In  the  collision  operator,  how¬ 
ever,  the  errors  are  exactly  the  same  from  step  to  step. 
In  addition  we  expect  that  the  radio  frequency  inhomo¬ 
geneity  leads  to  strongly  correlated  errors  in  the  lattice 
encoding.  Hence,  we  have  proposed  replacing  a  single 
collision  operator  with  a  set  of  collision  operators  that 
have  similar  fidelity  but  randomized  error  terms. 

Since  the  collision  operator  for  the  Burgers  equation  is 
a  zero-order  coherence  term,  the  collision  operator  com¬ 
mutes  with  the  rotation  operator.  Therefore,  we  apply 
a  90°  rotation  operator,  to  the  collision  operator  at  each 
step  to  mitigate  error  growth.  Consequently,  a  dramatic 
improvement  is  observed  as  shown  in  Figure  2.  On  a 
logarithmic  plot,  the  simulation  results  fit  a  line  with  a 
slope  of  3/4.  If  the  error  terms  in  the  collision  operators 
were  totally  randomized  and  hence  followed  a  Gaussian 
distribution,  the  best-fit  regression  line  should  have  had 
a  slope  of  1/2.  The  deviation  between  our  simulation 
data  and  the  ideal  Gaussian  case  indicates  residual  sys¬ 
tematic  error  in  the  collision  operator.  In  a  future  study, 
we  may  use  strongly  mudulated  pulses  to  randomize  the 
error  terms.  The  experimental  number  densities  are  over- 
plotted  in  Figure  3  with  the  exact  analytical  solutions. 
Eight  successive  time  steps  of  the  quantum  algorithm 


were  implemented  on  16  two-qubit  sites.  An  improve¬ 
ment  of  our  present  experimental  approach  using  collision 
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FIG.  3:  The  experimental  data  sire  plotted  together  with  the  ana¬ 
lytical  solutions  for  8  time  steps  on  a  lattice  of  16  parallel  two-qubit 
QIPs.  Viscosity:  \  ~~ .  Experimental  NMR  data  (dots)  versus 
analytical  solution  (curves).  Randomizing  the  error  terms  in  the 
collision  operator  has  improved  the  experimentsil  results  dramati¬ 
cally. 


operators  with  modulated  phases  is  observed.  The  agree¬ 
ment  of  the  data  to  the  analytical  solutions  is  encourag¬ 
ing  and  suggests  that  totally  randomizing  error  terms  in 
the  collision  operator  may  offer  further  improvement. 

NMR  quantum  simulations  has  provided  an  alternative 
way  to  study  the  NMR  spectroscopic  implementations. 
From  the  simulation,  we  find  the  major  error  sources 
are  due  to  imperfect  control  of  the  quantum  spin  system 
and  the  Fourier  approximation  associated  with  setting  its 
magnetization  profile.  Our  proposed  method  for  convert¬ 
ing  the  systematic  errors  into  random  errors  is  effective. 
The  improvement  we  achieve  relative  to  the  previous  ex¬ 
periment  is  encouraging,  and  it  demonstrates  the  possi¬ 
bility  of  using  the  same  technique  in  future  studies.  The 
closeness  of  the  numerical  data  to  the  exact  analytical 
results  for  the  nonlinear  Burgers  equation  further  proves 
the  practicality  of  implementing  the  QLG  algorithm  us¬ 
ing  a  spatial  NMR  technique.  In  addition,  although  the 
limitation  of  the  Fourier  approximation  is  not  dominant, 
the  problem  of  precisely  initializing  a  lattice  of  QIPs  still 
remains  an  open  issue. 

We  thank  M.  Pravia,  N.  Boulant,  H.  Cho  and  Y.  Liu  for 
valuable  discussion.  This  work  was  supported  by  the  Air 
Force  Office  of  Scientific  Research,  along  with  DARPA, 
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